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We derive a 1/c-expansion for the single-particle density matrix of a strongly interacting time- 
dependent one- dimensional Bose gas, described by the Lieb-Liniger model (c denotes the strength 
of the interaction). The formalism is derived by expanding Gaudin's Fermi-Bose mapping operator 
up to 1/c-terms. We derive an efficient numerical algorithm for calculating the density matrix for 
time-dependent states in the strong coupling limit, which evolve from a family of initial conditions in 
the absence of an external potential. We have applied the formalism to study contraction dynamics 
of a localized wave packet upon which a parabolic phase is imprinted initially. 



PACS numbers: 05.30.-d,03.75.Kk 



I. INTRODUCTION 



One of the most attractive many-body quantum sys- 
tems nowadays has been introduced by Lieb and Lin- 
iffer in their landmark paper more than forty years ago 
[l| . The system is composed of N identical (5-interacting 
bosons in one spatial dimension and is referred to as a 
Lieb-Liniger (LL) gas. They have presented an explicit 
form of the many-body wave function for a homogeneous 
gas with periodic boundary conditions Jj, including 
equations describing the ground state and the excitation 
spectrum. In the strongly interacting "impenetrable- 
core" regime [2], such a one dimensional (ID) system 
is referred to as the Tonks- Girardeau (TG) gas; exact so- 
lutions in this limit are obtained by Girardeau's Fermi- 
Bose mapping 0. Following Ref. Yang and Yang 
Q have eliminated a possible existence of phase tran- 
sitions in the LL system by proving analyticity of the 
partition function. After many recent experimental suc- 
cesses d, d, @, 0, B HI in realization of effectively one 
dimensional (ID) interacting gases, from the weak up to 
the strongly interacting TG regime Q , the LL model has 
attracted considerable attention of the physics commu- 
nity. There is a clear reason for this; nontrivial quantum 
many-body systems are notoriously oblique to a quan- 
titative analysis, and therefore possibility of an exact 
treatment in particular cases, together with experimental 
realization, is of great value. Moreover, exact solutions 
can be useful as a benchmark for approximate treatments 
aiming to describe a broader range of physical systems. 

Even though exact LL many-body wave functions can 




'Electronic address: ' rpezer@phy.hr| 



be constructed in some cases (e.g., stationary (ll. [lol [TTI. 

ITp . Il6| or time-dependent wave functions 
120|), the calculation of observables (correla- 
tion functions) from such solutions usuall y p oses a major 
difficulty in practice HI, HI, Hi, 111,11, |13, HI, ill, 
HOjUIIjIHI- Various methods have been employed to over- 
come this difficulty including, for example, the quantum 
inverse scattering method (e.g^see [2l|, |20|) and quan- 
tum Monte Carlo integration [2^ . A recent discussion of 
several exact methods for the calculation of correlation 
functions of a nonequilibrium ID Bose gas can be found 
in Ref. 32]. In the TG limit, the momentum distribution 
can be analytically studied for a rin g ge ometry, and also 
for harmonic confinement (e.g., see |33l. |3^). Numerical 
methods for the calculation of the reduced single-particle 
density matrix (RSPDM) can be performed efficiently 
for various TG states (ground state, excited and time- 
dependent states, see Ref. [11] for hard-core bosons on 
the lattice, and Ref. [11] for the continuous TG model 

Ultracold atoms in ID atomic wave guides enter the 
strongly interacting regime at low temperatures, in tight 
transverse confinement, and with strong effective interac- 
tions [13, m, HI] . The correlations functions of a LL gas 
can in this limit be calculated by using 1/c expansions 
(e.g., see [H, Hi, 113) from the TG (c oo) regime. 
These calculations in the strongly interacting limit ex- 
ploit the fact that a bosonic LL gas is dual to a fermionic 
system [40], such that weakly interacting fermions cor- 
respond to strongly interacting bosons and vice versa 
[Tol ITTj . A strongly interacting ID Bose gas was stud- 
ied in Ref. [12] by using perturbation theory for the dual 
fermionic system. In Ref. [13], the dynamic structure 
factor was calculated for zero and finite temperatures. 

Here we calculate the 1 /c correction for the RSPDM of 
a Lieb-Liniger gas, which adds upon a recently obtained 
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formula for the RSPDM of a TG gas The method 
is derived by using the 1/c term of the so-called Fermi- 
Bose (FB) mapping operator introduced by Gaudin p7l |. 
The FB operator method provides us with exact time- 
dependent solutions of a LL model [l^ in the ab- 
sence of external potentials and other boundary condi- 
tions; it was recently used to study free expansion of a 
LL gas 20]. We derive an efficient numerical algorithm 
for calculating the RSPDM for time-dependent states in 
the strong coupling limit, which evolve from a family of 
initial conditions in the absence of an external potential. 
We employ it to study the evolution of a many-body wave 
packet with a parabolic phase imprinted at t = 0, which 
corresponds to focusing with a lens in optics^ technique 
experimentally feasible in ID atomic gases 

This work complements the studies of nonequihbrium 
dynamics of interacting Bose gases which have been ad- 
dressed by use of the LL model away from 16 18 23, 33 



and in the TG limit 0, |M El H E lilMSi 
The phenomenological relevance of these studies is un- 
derlined by the fact that nonequihbrium dynamics is ac- 
cessible experimentally 0, ll] ■ 

The paper is organized as follows: Section |lT] gives a 
detailed account of the formalism, where we outline the 
procedure for the calculation of the RSPDM. In section 
mil we present the example of a Bose gas initially in the 
ground state within a harmonic trap, effectively in the 
TG regime; at i = a parabolic phase is applied on the 
initial state and the harmonic potential is turned off. The 
subsequent dynamics leads to focusing of the cloud and 
local increase of the 1/c correction. Details of the deriva- 
tion and relevant mathematical identities are collected in 
the Appendices. 



II. FORMALISM FOR THE TIME-DEPENDENT 
STRONGLY-INTERACTING LIEB-LINIGER GAS 



Schrodinger equation for a noninteracting Fermi gas, 



^ a2 I 



dt 



then the wave function 
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t) = ^jAf{c) Jl [sgn(a;; Xr)- 
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obeys Eq. H]) as pointed out in Ref. [TtI (see also 
Refs. [l9l,'20]). Here, 7V(c) is a normalization constant, 
and the differential operator in the square brackets de- 
notes the Fermi-Bose mapping operator [l9l | . The deriva- 
tives do not act on any of the sign functions (the sign 
functions can be avoided by working in only one sector of 
the configuration space, e.g., for a;i < a;2 . . . < xjv, e.g., 
see Refs. [iM [lO]). We assume that ^f{xi, . . . ,XN,t) 
is normalized to unity. Equation ([3]) can be reorganized 
in a finite power series with terms of order l/c™, where 
TO = 0, 1, ... , A''(iV— 1)/2. By using the Fermi-Bose map- 
ping operator in this form, one obtains a systematic ex- 
pansion of the exact many-body wave function ipLL in 
the inverse interaction strength 1 /c. 

The expansion of the wave function ([3]) in the inverse 
coupling strength 1/c is particularly useful in the strong 
coupling limit, as it allows an approximate calculation 
of one-body observables contained within the reduced 
single-particle density matrix (RSPDM). In the Tonks- 
Girardeau limit, where c = oo, a formula for efficient cal- 
culation of the RSPDM has recently been derived [s^ . 
Here we generalize that result to include the 1/c term 
in the expansion. By keeping only the 1/c terms in the 
expansion, Eq. ^ reduces to 



The Lieb-Liniger model describes a system of bosons 
interacting via pointlike interactions; the many-body 
Schrodinger equation for the LL gas of N such bosons 
reads yy] 



.dtp 



LL 



dt 



1=1 



LL 



dx. 



J2 ^cSix, 

l<i<j<N 



'Xj)lpLL- (1) 



Here, ?/'ll(xi, . . . , xjv, i) is the time-dependent bosonic 
wave function, c is the strength of the interaction. For 
now, we assume the absence of any external potential and 
boundary conditions. Under these circumstances, the 
time-dependent LL model ([T|) can be solved by employing 
the Fermi-Bose transformation [13, [3] ■ This method can 
be applied, e.g., to exactly study free expansion from an 
initially localized state [20| . If ipF (2^1 , ... , xn , t) is an an- 
tisymmetric (fermionic) wave function, which obeys the 



iPll{xi, ■ ■ - .XN^t) ~ V7V(c) J]^ sgn(x-,- 

l<i<j<Af 
dtpF dlpF^ 



- sgn(a::; - Xr) 



l<r<l<N 



dxi dXr 



(xi,...,Xk ,t) 



(4) 



The first term is simply the TG gas wave function 0], 
while the second term gives the 1/c correction to the TG 
wave function when the coupling constant c is finite. This 
expression is the starting point for all results that will be 
derived in this paper. The RSPDM is defined as 



PLL{x,y,t) =N Jdx2 ■ ■ ■dxNtpLL{x,X2, ■ . .,XN,t)* 

X-tpLLiy,X2,...,XN,t). (5) 

On inserting Eq. (jH) into Eq. (O) we obtain a formal 
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expression for the 0{l/c) correction of the RSPDM: 



Pll{x, y, t) = Af{c)I{X) i^*p(x, X, t)^F (y, X, t) 
1 

+ - 
c 



-M{c) ^ i{X) sgn(xz - Xr) 

l<r<l<N 



+ sgn(a;; - Xr) iP*f{x, X, t) 



My,x,t) (6) 



fdijF 


dijF\ 1 


\ dxi 


J (y^X,t) 



0(l/c2 



Here, X — {x2, ■ ■ ■ ,X]y) and the integral operator i{X) 
is defined as: 



i{X) ^nJI dxn sgn( 



X - Xn)sgn{y - Xn). (7) 



The first term on the right hand side of Eq. ^ is the 
TG gas RSPDM [3^. It can be proven (as is done in Ap- 
pendix |A| that only the partial derivatives with respect 
to the first coordinate xi in Eq. ^ give a nonvanishing 
contribution. After eliminating the vanishing terms from 
Eq. ^ we are left with 



Pll{x, y, t) = M{c)ptg{x, y, t) + -7V(c) 

c 



N 



^i{X) sg-!v{x~xi) 

1=2 

sgn{y - xi)iP*f{x,X, t) 



( di>i. 



V dxi 

f dtpF 
\ dxi 



'4>F{y,x,t)+ 



(x,X,t) 



iv,x.t) 



0(l/c2 



1, 



Ptg{x, y, t) + -[t^{x, y, t) + 7^{y, x, t)*] + 0{l/c^). 
c 

(8) 

In Eq. ([8]) we have implicitly defined the quantity 
■q{x,y,t): 



N 



v{x, y,t) H^) ^sMx - xi] 

1=2 

^{N -l)i{X)sgn{x-X2) 



dxi 

dipF 
dxi 



'iJjF{y,x,t) 



^Fiy,x,t). 



(x,X,t) 



(9) 



By using the calculation presented in Appendix |^ it fol- 
lows that 



dx[ri{x, x, t) + rj{x, x, t)*] — 0, 



but also lowers it in others such that the integral over 
the terms of order 1/c is zero. By using this result we 
see that in the leading 1/c approximation we can take 
M{c) — 1; this fact has already been utilized in the last 
line of Eq. ([8]). It is straightforward to verify that the 

integrals J(X)sgn(a; — xi) ( ^f- ) i'Fiy, X) are inde- 

pendent oil {2 < I < N), which yields the second identity 
in Eq. Q. The last line of Eq. ^ verifies that to or- 
der 1 /c the RSPDM possesses as required the symmetry 
PLL{x,y,t) = pLL{y,x,t)*. 

Up to this point we did not make any assumptions 
on the structure of ipF (except that it is antisymmet- 
ric and normalized), that is, the derivation was general, 
valid even if the wave functions should describe the gas 
in an external potential etc. Quite generally, ipp can also 
be considered as a function of c, and one could expand 
it in powers of 1/c. ^From this point on, ^f will be 
represented as a Slater determinant formed from single 
particle wave functions. 
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p 



{-)^(l)Pi{xi,t) ■ ■ ■(l)pNixN,t); (10) 



such wave functions can be used to study dynamics on an 
infinite line, which arises from initial conditions given by 
Eqs. (HI) and In Eq. (j}j{x,t) {j = 1,...,N) are 
orthonormal single-particle wave functions which obey 
the Schrodinger equation id(j)j/dt — —d^ifij/dx^. P de- 
notes a permutation P(l, . . . , N) = {PI, . . . , PN) of the 
particle number indices, and (— )^ is its signature. In 
this form ipF enables us to rewrite Eq. ^ such that it 
involves ID integrals only, resulting in certain algebraic 
cofactors suitable for numerical calculation. 



A. Algorithm for RSPDM calculation 

For the sake of clarity of the presentation, we first de- 
scribe the algorithm for calculation of the RSPDM, and 
only afterwards provide its derivation. Without loss of 
generality we consider the case x < y. The first step is 
to calculate the integrals 

Ikj{x,y,t) ^ 6ki ~2 dx (l)l{x' ,t)4>i{x' ,t) (11) 



and 



that is, the 1/c correction to the single-particle density 
Pll{x, X, t) increases the density in some regions of space, for fc, Z = 1 



hAy^t)^-5ki+2l dx'cj)lix',t)cj)iix',t); (12) 
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These integrals are arranged in the foUowing matrices: 



Fix, y,t) 



h,i{x,y,t) Ii^2{x,v,t) ... Ii,N{x,y,t) 
hAx,y,t) l2,2{x,y,t) ... l2.N{x,y,t) 

lN,i{x,y,t) lN,2[x,y,t) ... lN,N{x,y,t) 



and 



P«(a;,y,t) 



h.iix,y,t) hAx,y,t) ... Ii,i{y,t) ... Ii^N{x,y,t) 
l2,i{x, y,t) l2,2{x,y,t) ... I2,i{y,t) ... l2,N{x,y,t) 

lN,i{x,y,t) lN,2ix,y,t) ... lN.iiy,t) ... lN.N{x,y,t) 
I 



(13) 



(14) 



Let us define the column vector 

/ (j)i{x,ty 

*(x,t)= : 

\cj)N[x,t)^ 

and its first spatial derivative 

(c^{{x,ty 



(15) 



The TG reduced single-particle density matrix is given 
by M 



PTG{x,y,t) = det [P{x,y,t)] 

^\x,t)[Pix,y,t)~']'^^{y,t). (17) 



*'(x,t) 



(16) 



\(l)'^{x,t)^ 



The quantity r]{x,y,t) can also be written in a convenient matrix form (suitable for efficient numerical implemen- 
tation): 



TV 



1=1 



77(x,2;,i)-^ det P(')(^,y,t) *'t(x,t) P»(x,2;,i)-i ^{y,t) 



- det[Pix,y,t)]^'^ix,t)[Pix,y,t)-'] *(y,t) 
I 



(18) 



If any of the matrices P*^'^ happen to be singular, we can 
resort to a direct calculation via algebraic cofactors (see 
the proof of the algorithm in Appendix [B| . This hap- 
pens rarely and only for some particular high-symmetry 
points. Equations (fT3)) -(fT8 l) provide the grounds for an 
efficient numerical method for calculating the RSPDM, 
which is a generalization of the previously introduced 
method for the TG gas 36]. 

It is convenient to calculate the diagonal correction 
to the RSPDM. In this case, the matrices appearing in 
Eq. are very simple and it is straightforward to ob- 
tain 



PLL{x,X,t) 

Af{c) 



PTG{x,X,t) 

[Tr{p')Tr{l)-TT{p' -1)]^^ ,^ (19) 

where on the right hand side p is N x N matrix given by 

Pk,i{x,t) = (l)l{x,t)(l)i(x,t) 



and primes denote the first spatial derivative. We can 
also omit the normalisation constant from Eq. (|19p since 
A/'(c) = 1 within the 1 /c approximation considered here. 



III. CONTRACTION DYNAMICS OF 
STRONGLY INTERACTING ID BOSE GAS 

In this section we apply the presented formalism to 
investigate time-evolution of a many-body wave packet 
which contracts in space after a parabolic phase was im- 
printed on it. More specifically, we consider N strongly 
interacting bosons which arc initially in the ground state 
in the presence of a harmonic potential. We assume that 
this initial state is well described by the Tonks-Girardeau 
ground state wave function, a symmetrized Slater deter- 
minant. As described in more detail below, at time t = 0, 
a parabolic phase is suddenly imprinted onto this state, 
and the harmonic potential is turned off. Imprinting the 
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parabolic phase is experimentally feasible ^ and is equiv- 
alent to the action of a focusing lens on the travelling- 
wave state of a light beam in optics. In the context of 
atomic gases this can be achieved by applying, over a 
short period of time, a tight longitudinal harmonic po- 
tential to the one-dimensional gas. This results in pro- 
viding different initial momenta to different parts of the 
wave packet, with the momenta being proportional to the 
distance from the center of the wave packet and directed 
towards the center. 

In order to connect the model written in Eq. ([T]) to 
physical units we first write the ID Schrodinger equation 
for this system: 



dip 



LL 



dr 



N r 

E 

1=1 



■- - I HO c2 



LL 



E 

l<i<j<Ar 



where g is the ID scattering length [37[, ^ is the spa- 
tial coordinate, and t denotes time. Here we denote the 
axial trapping frequency by ojho- By multiplying the 
above expression with 2 / Hujho we obtain a dimensionless 
equation suitable for numerical calculation and a physical 
interpretation of the relevant scales: 



LL 



dt 



N 

E 



E 



■IpLL 

2.9 



l<i<i<JV 



aHoEno 



8{xi - Xj) 



LL- 



(20) 



where we have denoted zero point energy of the har- 

The unit of time is 



hijj 



HQ. 



monic trap with Eho = 
2/ujhoi while space is in units auo = ^Jf^/iTi^HO] in 
other words Xi — Ci/O'HO, and t = ojhot/2. If we 
set c — 2g / anoEHO we easily verify that it is exactly 
the interaction strength of the LL model c, see Eq. ([Ij. 
For concreteness, let us assume that the harmonic trap 
has frequency ijJho = 60 Hz, which is loaded with ®^Rb 
atoms, N = 12. These values yield ano = 1.38/im. As 
we have already stated, we assume that the effective in- 
teraction strength c is sufficiently large such that the sys- 
tem is initially in the Tonks-Girardeau regime, and that 
the 1/c correction is negligible a,t t — 0. Under these 
conditions the ground state is given by 



i^g.s. = 



n ^s'^^^j 

l<i<i<Af 



Xi) 



N 

det [(j)j{x,j 

j,m=l 



(21) 

where (f>j{x) is the jth eigenstate of the SP harmonic 



oscillator 



, ^ — j^j^j. At t = we imprint a 
parabolic phase to the system. Technically, each of the 
SP wave functions in the determinant (|2ip is multiplied 

by 



P{x) 



where we have chosen b = \f2. Thus, the initial wave 
function which starts to evolve in the absence of a trap- 
ping potential is 



i}LL{x\, ■ ■ ■,XN,0) 



n sgn(a;j 

l<i<j<N 



Xi) 



N 



det [(P,iXrn)P{Xrn)]/VNl (22) 
j,m=l 

After this initial excitation the wave packet starts to con- 
tract. Figure [1] displays the evolution of the SP x-space 
density and the momentum distribution (SP fc-space den- 
sity) for c = 50; the momentum distribution is defined as 



dx / dx' pLL{x,x',t). 



riLLik, t) = (27r)' 



It should be emphasized that the results are scalable, that 
is, the same results can be obtained for smaller values of 
c provided that the initial state is broader (this corre- 
sponds to the fact that for the LL gas in equilibrium 
the regime is determined by the interaction strength c 
divided by the linear density [H). As the gas becomes 
more dense during the contraction, the 1/c correction 
to the TG density becomes larger and clearly visible in 
the x-space density. Note that the gas compresses more 
strongly for the finite c in comparison to the TG results. 
After the wave packet reaches the maximally dense con- 
figuration, it starts to freely expand and the 1/c correc- 
tion starts to become less visible. It is interesting to 
note that while the changes in the x-space density are 
nonnegligible in contraction, the 1/c correction in the fc- 
space density is practically negligible at all times. ^From 
our simulations we find quite generally that the 1/c cor- 
rection to the fc-space density is much less sensitive to 
changes in the coupling strength than the x-space density. 
Note that the total x-space density in the presented form 
of the 1/c approximation can become negative in some 
parts of the space, if the approximation is used beyond 
its range of validity (e.g., if c is not sufficiently large). 
This gives a clear indication that higher order terms, of 
order 0{l/c^), become relevant, at least in those parts 
of the space where the density is negative. 

We emphasize again that the results presented here, 
for an evolution in one dimension without further bound- 
ary conditions, are correct up to 1/c for all times of the 
evolution, i.e., there is no accumulation of the error in 
time due to the approximation and we can extract the 
asymptotic behaviour of the LL gas. 



IV. CONCLUSION 

We have studied the 1/c expansion in the strong cou- 
pling limit of the Lieb-Liniger model, exploiting the for- 
malism based on a mapping between free Fermi and inter- 
acting Bose gas. More specifically we have developed an 
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FIG. 1: (color online) The a;-space density (left, red solid line), and momentum distribution (right, red solid line) of the LL 
gas with N=12 bosons and c = 50 at different times of the propagation. At the time of 1.3 ms the cloud attains the maximal 
compression. The corresponding shapes for the TG gas are also shown as black dot-dashed lines. The difference between 
red solid and black dot-dashed lines corresponds to the 1/c correction. The curves are shifted on the ordinate axis for better 
visibility. See text for details. 



efficient algoritlim to calculate tlie reduced one body den- 
sity matrix of tlie strongly interacting LL gas, which is 
based on the metfiod previously developed for the Tonks- 
Girardeau gas ^3^, and the Fermi-Bose transformation 
[17, 19, 20]. The method is suitable to describe dynamics 
on an infinite line, in the absence of external potentials, 
and for initial states which can be described by Equations 
dH) and pi7)) . An attractive feature is stable accuracy for 
all times since the error is not being accumulated dur- 
ing the evolution. For the x-space density we provide 
a very compact formula Eq. (|19[) that can be employed 
for quite large numbers of particles enabling fast calcu- 
lation of the nonequilibrium x-space density dynamics. 
In addition, the efficiency of the method is demonstrated 
in an example discussed in Sec. IIIII we have examined 
nonequilibrium dynamics where a parabolic phase is im- 
printed onto a localized wave packet after which the wave 
packet contracts for a while. Our simulation studies the 
setup which was experimentally realized in Ref. Q for a 
weakly interacting Bose gas at finite temperature. 



tion for Science. 



APPENDIX A: FORMULAS 

Integrals in Eq. ^ involving derivatives in integration 
coordinates all vanish due to following expression for k ^ 
1 and for any I ^ k: 



I{X)iig\i{xk - xi) 



\ dxk J 



ix,X) 



+ i{X)agn{xk - xi)ip*p{x, X) 



d 



dxk 



= /(X)sgn(xfc - a^O^ [rF{x,X)^F{y,X)] 

This is obtained by collecting the corresponding terms 
from the first and the second part of the 1/c correction 
to the RSPDM in Eq. (©. Let us now look at the fcth 
integration only. Integrating by parts we have: 
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[sgn(a;fc - xi)il)*p{x, X)'il)F{y, X)] 

c 

dxk 2S{xk - xi)il)*p{x, X)'ipF{y, X) 



The first term is zero due to boundary condition at infin- 
ity and the second term is zero because of the antisym- 
metry of tAf- By using the method presented above it is 
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straightforward to see that 

dx [ri{x, X, t) + ri{x, x, t)*] = 



where 77 is defined in Eq. ([9|). A direct consequence 
is that the normahzation constant A/'(c) = 1 in the 1/c 
approximation cosidered here. 



APPENDIX B: DETERMINANT EXPANSION 



In this appendix we derive Eq. (|18|) . The derivation 
is similar to that for the RSPDM of a TG gas To 
simplify the notation, we will suppress the time variable 
in all formulas. By inserting Eq. (jlOp for ^/^j? in Eq. ^ 
for rj we obtain 



(Bl) 



Here, P G S^-i is a permutation of the — 1 indices 
not including i: (p2, . . • ,Pn) = -P(l ...i— 1 i + 1 . . . N). 
Similarly, Q G Sn-i permutes the N — I indices not 
including j: (52, . . . ,qN) = Q{1 ■ ■ ■ j - I j + I ■ ■ ■ N). 

We consider the sum over the permutations. 
Let {q'^ , . . . , q'^) be the ordered series of integers 
{1, . . . ,N}\{l,j}, where / ^ j, and Q' G Sn-2 a 
permutation of these N — 2 numbers, ((73,..., g^) = 



Q'{q'^, . . . , q'pf)- With these, the sum over permutations 
in Eq. (|B1[) can be written as 



1 \"(-)P+Q7 I .../ 



{N-2) 



N 



f7V-2)l^^ ) ^P2-J^P3,q'3- ■ ■ 1} 

^ ''1=1 P,Q' 



N 



PiV,<?N 



( = 1 P 



J2 detP\[>{x,y) 



(B2) 



Here, j isaiV— IxiV— 1 matrix obtained by elim- 
inating the ith row and the jth column from the matrix 
P*^'^ defined in Eg. p^ . This yields the following result 
for ri{x, y): 

V{x,y) ^Y.^^-{x)r^,{y) ^ (-r+^ detPg(x,j/) 

ij 1 = 1 

(B3) 

The final result, Eq. (fT5|) . is obtained from (jB3p with the 
use of standard matrix algebra connecting minors, the 
cofactor matrix, and the matrix inverse. 
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